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Abstract. The Kluzniak & Abramowicz model explains high frequency, double peak, "3:2" QPOs observed in 
neutron star and black hole sources in terms of a non-linear parametric resonance between radial and vertical 
epicyclic oscillati ons of an almost Kepl erian accretion disk. The 3 : 2 ratio of epicyclic frequencies occurs only 
in strong gravity. lRebuscc-1 <l2004l) and iHorak I (|2004) studied the model analytically: they proved that a small 
forcing may indeed excite the parametric 3:2 resonance, but they have not explained the physical nature of the 
forcing. Here we integrate their equations numerically, dropping the ad hoc forcing, and adding instead a stochastic 
term to mimic the action of the very complex processes that occur in disks as, for example, MRI turbulence. We 
demonstrate that the presence of the stochastic term triggers the resonance in epicyclic oscillations of nearly 
Keplerian disks, and influences their pattern. 

Key words. Methods: data analysis - Methods: statistical-X-ray: binaries-Relativity-Accretion, accretion disks 



'O ! 1. Introduction l2003tlLi &: Naravan ||200^ . The Kluzniak & Abramowicz 

j_i ■ resonance model ( see a collection of review articles in 

03 ; Quasi Periodic Oscillations (QPOs) are a common phe- lAbramowiczl 120051) stresses the importance of the ob- 

nomenon in nature. In the last few years many kHz QPOs served 3:2 rat i , pointing out that the commensurability 

have been detected in the light curves of about 20 neu- of frequencies is a clear signature of a resonance. The rel- 

tron star and few black hole sources (for a recent review, evance of the 3 : 2 ratio and its intimate bond with the 

see |van der Klis | |200J). The nature of these QPOs is q POs fundamental nature is supported also by recent ob- 

one of the mysteries which still puzzle and intrigue astro- ser vations: Jeroen Homan of MIT reported at the AAS 

physicists: apart from giving important insights into the mee ting on the 9th of January 2006 that the black hole 

disk structure and the mass and spin of the central object candidate GRO J1655-40 showed in 2005 th e same QPOs 

(e.g. lAbrarriowicz fc Kluzniak I |200lj |Aschenbach | |200J (at 300 Hz and ^ 450 Hz) first detected bv lStrohmaverl 

iTorok et al. Il2005|) . they offer an unprecedented chance J200ll ). 
to test Einstein's theory of General Relativity in strong 

fields. The main limitation of the resonance model is that 

High frequency QPOs lie in the range of orbital it does not yet explain the nature of the physical mech- 

frequencies of geodesies just few Schwarzschild radii anism that excites the resonance. The idea that turbu- 

outside the central source. This feature inspired sev- lence excites the re sonance and feeds energy into it (e.g. 

eral models based directly on orbital m otion (e.g. lAbramowicz 1 12005|) is the most natural one, but it has 

IStella fc Vietril Il998l lLamb fc MilleTI l2003lh but there never been explored in detail. The turbulence in accre- 

are als o models that are based on accretion disk oscil- tion disks is most p robably due to the Magne to-Rotational 

lations ijWagoner et al. ll200lilKatc~ll200lilRezzolla et al. I Instability (MRI; iBalbus &: Hawlev"1ll99lj) . At present, 
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numerical simulations of turbulence in accretion disks do 
not fully control all the physics near the central source. 
For this reason, they cannot yet address the question of 
whether MRI turbulence does play a role in exciting and 
feeding the 3 : 2 parametric resonance. A situation like 
this is not specific of astronomy, but it is shared by other 
fields in applied research and engineering. The most com- 
mon and, at the same time, effective, solution consists of 
modelling the unknown processes as stochastic ones. Such 
processes are characterized by a huge number of degrees 
of freedom and therefo re they can be assumed to have a 
stochastic nature (e.g. iGarcia- Okialvo fc Sancholll999h . 
Lacking any a priori knowledge, the most natural choice is 
represented by Gaussian white-noise processes. Of course, 
such an assumption is only an approximation. However, 
it can provide an idea of the consequences on the sys- 
tem of interest of the action of a large number of complex 
processes. This approach leads to the modelling of phys- 
ical sys tems by means of stochastic differential equations 
(SDE) JMavheck lITflTfll Xv^k lOhauem fc Snauos I ITflflTl 
lOa.rcia-Okialvo fc Saucho UToOflt IVio et al. 11200^. 

The present paper is a first qualitative step in this 
direction in the context of QPO modelling. In Sec. 
we synthesize a stochastic version of the non-linear res- 
onance model. Some experiments are presented and dis- 
cussed in Sec.|3j The last section summarizes our findings. 
Since SDEs are not yet very well known in astronomy, 
Appendix lAl provides a brief description of the techniques 
for the numerical integration that are relevant for practical 
applications. 

In all the experiments, we adopt the units tq = 
2GM/c 2 = 1 and c = 1. 



2.2. Dynamics of a test particle 

A simple ma thematica l appr oach to this i dea w as first 
developed bv lRebuscol l|2004 and lHorakl l|2004 - in the 

context of isolated test particle dynamics. 

The time evolution of perturbed nearly Keplerian 
geodesies is given by 

z(t)+u%z(t) = f\p(t),z(t),ro,6 ]-, (1) 
p(t) + Lo 2 rP {t) = g[p(t), z(t), r , 6o]- (2) 

Here p(t) and z(t) denote small deviations from the cir- 
cular orbit r( h 0Q (radial and the vertical coordinates re- 
spectively), / and g account for the coupling, and cug 
and uj r are the epicyclic frequencies. In the case of the 
Schwarzschild metric, a Taylor expansion to third order 
leads to: 

f(p,z,r o ,0 o ) 

= cnzp + c b zp + c 2 ip 2 z + ci b pip + c 03 z 3 ; (3) 
g(p, Z, r Q ,9 ) = e Q2 z 2 + e 20 p 2 + e z2 z 2 + e 3Q p 3 

+ ei ze2 pz 2 + e 12 pz 2 + e r2 p 2 + e lre2 p 2 p. (4) 

The functional f orm of the coefficients Cj and ej can 
be found in iRebusco I l|2004j) . They are constants, which 
depend on ro, the distance of the unperturbed orbit 
from the centre. In previous studies these non-linear 
diffe rential equations hav e been integrated numerically 
llAbramowicz et al. I EoO^ and analyzed through a per- 
turbative method. These coupled harmonic oscillators dis- 
play internal non-linear resonance, the strongest one oc- 
curs when u)g : lo t — 3 : 2 and the observed frequencies 
are close (but not equal) to the epicyclic ones. 



2. A simplified model for kHz QPOs 

2.1. The Kluzniak-Abramowicz idea 

The key point of th e mec hanism proposed by 
lAbramowicz fc Kluzniak is the observation 

that kHz QPOs often occur in pairs, and that the 
centroid fre quencies of these pairs are in a rational 
ratio (e.g., IStrohmaver I l200l|) . This feature suggested 
to them that high frequency QPOs are a phenomenon 
due to non-linear resonance. The analogy of radial and 
vertical fluctuations in a Shakura-Sunyaev disk with the 
Mathieu equation pointed out that the smallest (and 
hence strongest) possible resonance is the 3 : 2. In all 
four micro-quasars which exhibit double peaks, the ratio 
of the two frequencies is 3 : 2, as well as in many neutron 
star sources. Moreover, combinations of frequencies and 
sub-harmonics have been detected: these are all signatures 
of non-linear resonance. A confirmation of the fact that 
kHz QPOs are due to orbital oscillations comes from the 
scaling of the frequen cies with 1/M, where M is the m ass 
of the central object l|Mc01intock fc Remillard"ll2004l) . 



2.3. Additional terms 

As we have seen the perturbation of geodesies opens 
up the possibility of internal resonances. However these 
epicyclic oscillations would not be detectable without 
an y source of energy to ma ke t heir ampli t udes grow. 
In lAbramowicz et al. I l)2003|) and IRebusco I l|2004|) this 
source of energy was inserted by introducing a parame- 
ter a. The effect of forcing (e.g., due to the neutron star 
spin), and its potential to produce new (external) reso- 
nances, have been addressed recently (e.g. Abramowicz 
2005). The main l imit in t he approa c h pro posed by 
lAbramowicz et al. I ((2003) and IRebusco I l|2004j) is that it 
represents an ad hoc solution. Moreover, as stressed in 
Sec. H it does not consider the many processes that take 
place in the central region of a n accretion disk as, for ex - 
ample, MRI-driven turbulence ijBalbus fc Hawlev1ll99l|) . 
For this reason, we propose the stochasticized version of 
Eqs. (UJ-© 

2(f) + wf *(t) - f[p(t),z(t),r , 6 ] = a z p; (5) 
p{t)+Lu 2 p(t)-g[p(t),z{t),r ,6o} = 0, (6) 

with o~ z a constant and /3(f) a continuous, zero mean, unit 
variance, Gaussian white-noise process. 
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Table 1. Resonant frequencies (apart from the epicyclic 
ones) for different initial conditions (1,2,3) and noise stan- 
dard deviation (er 2 = 0, 1CT 5 , 1CT 4 ). 



There is no full understanding of turbulence in accre- 
tion disks. We know that the radial component is fun- 
damental in producing the effective viscosity which al- 
lows accretion to occur, and that MRI-turbulence should 
be different in the vertical and radial direction. Here 
we make a first step by introducing a noise term only 
along the vertical direction: in the end this ansatz alone 
gives interesting results. In t he Shakura & Sunyaev model 
l|Shakura fe Sunvaev I Il973|) the turbulent viscosity is 
parametrized via the famous ass- It is reasonable to as- 
sume that a z is at maximum a fraction, smaller than ass> 
of the disk height. Hence for a geometrically thin disk one 
would expect a maximum a z ~ Kr 4 -1CT 3 . As shown in 
Appendix lAl the smallness of the stochastic perturbation 
permits the development of efficient integration schemes 
for the numerical integration of the system J^}-©. 

3. Results 

We explored the dynamics of the test particle for different 
values of a z and initial conditions z(0) and p(0). All the 
integrations are performed by means of the scheme Q A. 14)1 - 
with h = 5 x 1(T 4 and t = 10 5 , for r = 27/5 which 
is the value for which the unperturbed frequencies are in 
a 3 : 2 ratio. As a sample, three different starting val- 
ues [z(0),i(0),p(0),/j(0)] have been used: [0.01,0,0.01,0], 
[0.1, 0, 0.1, 0], and [0.2, 0, 0.2, 0], which we refer to as mod- 
els 1, 2 and 3 respectively. For each of them we considered 
three values of a z : 0, 10~ 5 and 10 -4 . As pointed out in the 
previous section, a noise level stronger than these is un- 
like to occur, since it would destabilize the accretion flow. 
Moreover when the initial perturbations z(0) and p(0) are 
greater than about 0.5 they diverge, even in absence of 
noise: this is the limit for which the system can be con- 
sidered weakly non-linear and physically meaningful. 

The lower panels of Figs.n-EDshow how the amplitudes 
reach greater values for greater noise dispersion. These 
plots are done for the initial conditions 2, but similar be- 
havior is also obtained for different initial conditions: as 
expected, noise triggers the resonances. With regard to 
the frequencies at which the resonances are excited, the 
dominant one are always the epicyclic frequencies (the 
strongest peaks in the upper part of the plots). However, 
the sub- and super-harmonics also react (see Tab. 1), and 
their signal is stronger for greater noise dispersion. As pre- 
dicted by means of the perturbative method of multiple 
scales, the dominant oscillations have frequencies (u* and 
uig), close to the epicyclic ones. The pattern of the other 



resonances (Tab. 1) is not interesting in itself, as it de- 
pends on the initial conditions and on the noise, but it 
is significant from a qualitative point of view, as it is a 
signature of the non-linear nature of the system. 

When the noise is ~ 10~ 3 or greater the solution di- 
verges, whilst when it is too small (~ 10 -6 ) it does not 
differ too much from the results without noise. The exact 
limit of o z over which the epicycles are swamped depends 
on the initial conditions: it is indeed lower for greater ini- 
tial conditions, and vice versa. 

In the case where noise is assumed to be due to MRI 
turbulence, this simple experiment constrains its ampli- 
tude: turbulence that is too low does not supply enough 
energy to the growing resonant modes, whilst too much 
turbulence prevents the quasi-periodic behavior from oc- 
curring. From this oversimplified model we get an indica- 
tion that the standard deviation of vertical MRI must be 
~ 10~ 5 — 10~ 4 , which is reasonable since it is comparable 
with a small fraction of the disk height. 

In a yet u npublished work (private communication, 
ISkinner1l2005j) considers how far the data from a QPO 
source can constrain the properties of a simple damped 
harmonic oscillator model - not only its resonant frequency 
and damping but also to some extent the excitation. Not 
unlike the present work, he adds random delta function 
shots to a simple harmonic oscillator equation, changing 
the amplitude and frequency of shots. He observes that 
the data constrain the allowed range of parameters for 
the excitation. 

4. Conclusions 

Up to now models for kHz QPOs have been based on de- 
terministic differential equations. The main limits of these 
models is that they correspond to unrealistic physical sce- 
narios where the many and complex processes that take 
place in the central regions of an accretion disk are not 
taken into account. In this paper, we have partially over- 
come this problem by adopting an approach based on 
stochastic differential equations. The assumption is that 
the above mentioned processes are characterized by a huge 
number of degrees of freedom, hence they can be assumed 
to have a stochastic nature. In particular, we have inves- 
tigated a simplified model for the Kluzniak-Abramowicz 
non-linear theory and shown that a small amount of noise 
in the vertical direction can trigger coupled epicyclic oscil- 
lations. On the other hand too much noise would disrupt 
the quasi-periodic motion. This is si milar to the stochas- 
ticall y excited p-modes in the Sun l|Goldreich fc Keelev I 
Il977h . 

From our simple example we get an indication that 
the standard deviation of vertical noise cannot be greater 
than 10~ 5 - 10~ 4 r„, nor smaller than ~ 10 6 r g , but 
better modelling needs to be done. Nonetheless good es- 
timates are still possible without detailed knowledge of 
all the mechanisms in accretion disks; this approach has 
the power to lead to a better understanding of both kHz 
QPOs and other astrophysical phenomena. 
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Fig. 1. Numerical simulation of the system J^J-JSJ. The 
upper panels show the power spectra of z(t) and p(t), 
whereas the lower panels show the corresponding phase 
diagrams. Here a z = (i.e. noise- free system). The dis- 
placements are in units of r g , the frequencies are scaled to 
kHz (e.g. assuming a central mass M of 2Mq). 

Appendix A: Some notes on the numerical 
integration of SDEs 

A.l. General remarks 

A generic system of SDEs can be written in the form 

x = a(t,x) + i:(t,x)P (A.l) 

Here x, a are n-dimensional column vectors, (3 is a Tri- 
dimensional column vector containing zero mean, unit 
variance, Gaussian white-noise processes and £ is a n x m 
matrix. Typically, this equation is written in the more rig- 
orous form 

dx = a(t,x)dt + H(t,x)dw, (A. 2) 

with solution 

x t =x to + / a(s,x s )ds + / H(s,x s )dw s . (A.3) 
Jt Jto 

Here to is a m-dimensional Wiener process. The numerical 
integration of SDEs is quite a difficult problem. In fact, in 



Fig. 2. The same as in Fig. G] but with er 2 = 10 5 . 



the case of ordinary differential equations (ODEs) 

dx = a(t, x)dt (A. 4) 

numerical integration schemes are, either directly or indi- 
rectly, based on a Taylor expansion of the solution 

x t = x tg + / a(s,x s )ds. (A. 5) 

Jto 

Something similar holds also for SDEs. However, the 
stochastic counterpart of the deterministic Taylor ex- 
pansion is rather complex. In order to understand this 
point without entering into overly technical arguments, 
it is instructive to compare the expansions relative to a 
one-dimensional autonomous version of Eq. l|A.4fl and of 
Eq. HA. 2D with m = 1. In this case, for the ODE l|A.4jl 
the first-order integral form of the Taylor formula in the 
interval [to, t] is 

xt = xt + a(x to ) / ds + R 2 , (A. 6) 

Jt 
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Fig. 3. The same as in Fig. Q] but a z — 10~ 4 . The com- 
parison with the phase diagrams in the previous plots in- 
dicates how much the amplitudes grow under the effect of 
slightly stronger noise. 

where R 2 is the remainder. For the SDE i|A.2|) the corre- 
sponding expansion is 

xt — xt + a i x to) / ds + a(x to ) / dw s (A. 7) 

+ a' (x to )a(x to ) / / dw z dw s + R, (A. 8) 
Jto Jto 

where the symbol " ' " denotes differentiation with respect 
to x, and R is the remainder. From Eq. I|A.7|I it is possible 
to see the presence of the additional terms 



(A.9) 



to Jto 



When n,m 7^ 1, it is possible to show that in the higher 
order expansions some quantities appear as 



dwl 1 • ■ • dwi l ~] dwl' , 



'to JtQ Jto 

(A.10) 

where ji,j2,"',ji S [0,1, ... ,m]. Such quantities are 
termed multiple stochastic integrals. The main problem 



in dealing with them is that they cannot be computed ex- 
actly. Unfortunately, in its turn, numerical approximation 
is also a difficult affair. 

The consequence of this situation is that, even in the 
case of simple systems, only integration schemes of very 
low order strong convergence 1 can be used. In fact, for the 
autonomous version of system (|A.2fl the most commonly 
used technique is the Euler scheme 

(A.ll) 



-<[k] 



o.[k]h[k] 



£[ fe ]Au>[k], 



where = t[k+i] — i[fc] is the integration time step at 
the time in-i , and the elements of the vector 

rtk+i 



dw t = w t 



k + l] 



(A.12) 



are independent identically-distributed Gaussian random 
variables with mean equal to zero and variance equal to 
h [k] . 

A. 2. Small noise approximation 

If one takes into account that the order of strong conver- 
gence for the scheme (|A.12|I is only 7 = 0.5, in contrast 
to 7 = 1 for its deterministic counterpart, then it easy 
to understand why SDEs are not yet a standard tool in 
physical applications. 

In order to im prove this situation, 
iMilstein fc Tretvakov I ljl997|) note that in many problems 
the random fluctuations that affect a physical system are 
small. This means that the system l|A.2(l can be written 
as 

dx = a(t,x)dt + eE(t,x)dw, (A. 13) 

where e is a small positive parameter. This is an impor- 
tant observation since, for small noise, it is possible to 
construct special numerical methods that are more effec- 
tive and easier to implement than in the general case. 
In fact, the term of the expansion depends not only on 
the time step h but also on the parameter e. Typically, 
the mean-square globa l error of the schemes proposed by 
IMilstein fc Tretvakov I l|l997j) is of order O ( h? + e k hi ) with 
< q < p. Although the strong order of these methods is 
given by q, typically not a large number, they are able to 
reach high exactness because of the factor e fc at h q . For 
example, the simple scheme 

x [k+ i] = x [k] + -{K x + 2K 2 + 2K 3 + K 4 ) + e£A™ [fc] 

(A.14) 

where 



Ki = ha(t [k] ,X[ k ]), 
K 2 = ha(t [k] +h/2,x [k] +Ki/2), 
K 3 = ha(t [k] +h/2,x [k] +K 2 /2), 
K 4 = ha(t [k+1] ,x [k] +K 3 ), 



(A.15) 
(A.16) 
(A.17) 
(A.18) 



We shall say that a discrete time approximation irw con- 
verges strongly with order 7 > at time T if there exists a 
positive constant C, which does not depend on 5, and a <5o > 
such that e(<5) = E(\xt — X[ T ] |) < CP for each 5 G (0, So). 
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is of order 0(h 4 + eh + e 2 h}/ 2 ). In other words, the order 
of strong convergence is 0.5, as for the Euler scheme, but 
better results are to be expected because of the term e 2 
that multiplies h 1 / 2 . 
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